function maskind=autodetectroi(image,N,thres)
% maskind=autodetectroi(image)
% maskind=autodetectroi(image,N)
% maskind=autodetectroi(image,N,thres)
%
% Detects multiple regions of interest (ROIs) from given image.
% N is a parameter used for median and mean filters (NxN-mask, see
% medfilt2 and filter2 for details. Basically, larger N, more 
% smoothness) and thres is the threshold level of the noise.
% Default value for N is 5. The found masks are outputted as a cell.
% By default, the threshold level is chosen automatically for each
% of the found peaks. If given, the threshold is applied globally.
%
% 2012 Ari-Pekka Honkanen
%

	Neighbors=10;	%Maximum radius to look local maximas, also defines edge width

	try,	N;
	catch,	N=5;
	end
	try, manual;
	catch, manual=0;
	end

	%Removes the data at the edges of the detector 
	%(the local maxima algorithm can't reach there, so this removes troubles experienced by the magic wand tool) 
	img=image;
	img([1:N end-N:end],:)=0;
	img(:,[1:N end-N:end])=0;

	%Apply a median filter to the data
	img=medfilt2(img,[N,N]);	

	%Apply a mean filters the data (smoothes edges)
	h=fspecial('average',[N N]);
	img=filter2(h,img);	

	try,	isthres=1;thres;
	catch,	thres=autogetthreshold(img);isthres=0;display(sprintf('Automatic threshold at %f.',thres));
	end

	imgthres=img;
	%Threshold the data
	for ii=1:size(imgthres,2)		
		imgthres(find(imgthres(:,ii)<thres),ii)=0;
	end

	%Finds the local maxima
	[peaks(:,1) peaks(:,2)]=imagemaximas(imgthres,10);

	%Remove maxima that are too close to each
	ii=1;
	while(ii<size(peaks,1)+1)
		for jj=size(peaks,1):-1:1
			if(jj~=ii)
				if((peaks(ii,1)-peaks(jj,1)).^2+(peaks(ii,2)-peaks(jj,2)).^2<10.^2)
					peaks(jj,:)=[];
				end
			end
		end
		ii=ii+1;	
	end

	%Create a raw mask (for manual threshold)
	if(isthres)
		img=imgthres;
		for ii=1:size(img,2)	
			img(find(img(:,ii)>=thres),ii)=1;
		end
	end

	%Apply teh magick!
	for ii=1:size(peaks,1);
		if(isthres)
			masks{ii}=magicwand(img,peaks(ii,1),peaks(ii,2));
		else
			masks{ii}=magicwand(img,peaks(ii,1),peaks(ii,2),img(peaks(ii,1),peaks(ii,2))/2);
		end
		for jj=1:size(masks{ii},2)
			masks{ii}(find(masks{ii}(:,jj)<0),jj)=0;
		end
	end

	%Remove masks that have overlap (in case of the overlap the masks should be the same)
	ii=1;
	while(ii<length(masks)+1)
		for jj=length(masks):-1:1
			if(jj~=ii)
				if(sum(sum(masks{jj}.*masks{ii}))~=0)
					masks(jj)=[];
				end
			end
		end
		ii=ii+1;
	end
	
	if 1	%plot image of found masks
		try
			masksum=zeros(size(masks{1}));
	
			for ii=1:size(masks,2)
				masksum=masksum+ii*masks{ii};
			end
			imagesc(masksum);
			title('The found masks')
		catch,end;
	
	end

	if(size(masks,2)==0)
		display('No masks found!');
	elseif (size(masks,2)==1)
		display('Found 1 mask!');
	else
		display(sprintf('Found %d masks!',size(masks,2)));	
	end

	%Converts masks to the new format
	for ii=1:size(masks,2)
		maskind{ii}=find(masks{ii});
	end	

end

function [row col]=imagemaximas(image,radius)

	N=1;

	for ii=1+radius:size(image,1)-radius
		for jj=1+radius:size(image,2)-radius
			if(image(ii,jj)>0&&image(ii,jj)==max(max(image(ii-radius:ii+radius,jj-radius:jj+radius))))
				row(N,1)=ii;
				col(N,1)=jj;
				N=N+1;
			end
		end
	end
end

function mask=magicwand(image,row,col,thres)

	try, thres;
	catch, thres=1;
	end

	mask=zeros(size(image));
	mask(row,col)=1;
	newpix=[row col];

	while(size(newpix,1)>0)
		aux=[];
		N=1;
		for ii=1:size(newpix,1)
			if(newpix(ii,1)>1&&mask(newpix(ii,1)-1,newpix(ii,2))==0) %Pixel not checked yet
				if(image(newpix(ii,1)-1,newpix(ii,2))>=thres)
					mask(newpix(ii,1)-1,newpix(ii,2))=1;
					aux(N,:)=[newpix(ii,1)-1 newpix(ii,2)];
					N=N+1;
				else
					mask(newpix(ii,1)-1,newpix(ii,2))=1;	%Enlargens the mask a little bit
				end
			end
			if(newpix(ii,1)<size(image,1)&&mask(newpix(ii,1)+1,newpix(ii,2))==0)
				if(image(newpix(ii,1)+1,newpix(ii,2))>=thres)
					mask(newpix(ii,1)+1,newpix(ii,2))=1;
					aux(N,:)=[newpix(ii,1)+1 newpix(ii,2)];
					N=N+1;
				else
					mask(newpix(ii,1)+1,newpix(ii,2))=1;
				end
			end
			if(newpix(ii,2)>1&&mask(newpix(ii,1),newpix(ii,2)-1)==0)
				if(image(newpix(ii,1),newpix(ii,2)-1)>=thres)
					mask(newpix(ii,1),newpix(ii,2)-1)=1;
					aux(N,:)=[newpix(ii,1) newpix(ii,2)-1];
					N=N+1;
				else
					mask(newpix(ii,1),newpix(ii,2)-1)=1;
				end
			end
			if(newpix(ii,2)<size(image,2)&&mask(newpix(ii,1),newpix(ii,2)+1)==0)
				if(image(newpix(ii,1),newpix(ii,2)+1)>=thres)
					mask(newpix(ii,1),newpix(ii,2)+1)=1;
					aux(N,:)=[newpix(ii,1) newpix(ii,2)+1];
					N=N+1;
				else
					mask(newpix(ii,1),newpix(ii,2)+1)=1;
				end
			end
			
		end

		newpix=aux;
	end

end
